Przykładowy przebieg doboru hiperparametrów

mammography_x = read.csv(file = 'mammography_x.csv', header=FALSE)
mammography_y = read.csv(file = 'mammography_y.csv', header=FALSE)

mammography = mammography_x %>% mutate(class = mammography_y$V1)

Aby lepiej zrozumiec charakterystyke danych, wyswietlmy je i zaznaczmy znane anomalie zanim przejdziemy do nienadzorowanej detekcji anomalii.

plot_ly(mammography,
        x = ~V1,
        y = ~V2,
        z = ~V3,
        color=~class,
        colors=c('#ffff00', '#000000'),
        opacity=0.7,
        size=1
      )

## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter3d

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

W przypadku danych o więcej niż 3 wymiarach, ciężko jest zwizualizować dane tak by widać było anomalie na pierwszy rzut oka.

Dowiedzmy się przynajmniej ile takich anomalii w zbiorze jest, aby poznać skalę problemu.

sum(mammography_y)

## [1] 260

Podział zbioru danych na treningowe i testowe z informacją o ilości anomalii w podzbiorach.

split_ratio = 0.8

set.seed(123)
split = sample.split(mammography, SplitRatio = split_ratio)
data_train = subset(mammography, split == TRUE)
data_train_x = data_train %>% select(1:6)
data_train_y = data_train %>% select(7)

data_test = subset(mammography, split == FALSE)
data_test_x = data_test %>% select(1:6)
data_test_y = data_test %>% select(7)

sum(data_train_y)

## [1] 185

sum(data_test_y)

## [1] 75

Grupowanie

K-Means

Do pogrupowania danych z użyciem algorytmu centroidów, trzeba wiedzieć ile takich centroidów ma być w danych. Sprawdzimy to z wykorzystaniem funkcji fviz_nbclust z pakietu factoextra. Jako, że można użyć kilku metod do sprawdzenia optymalnej ilości klastrów, sprawdźmy więcej niż jedną.

fviz_nbclust(data_train_x, kmeans, method = "silhouette")

fviz_nbclust(data_train_x, kmeans, method = "wss")

fviz_nbclust(data_train_x, kmeans, method = "gap_stat", k.max=7, nboot=30)

Zwizualizujmy, jak wygląda grupowanie dla wybranej ilości centroidów.

no_clusters = 3
cluster_data = kmeans(data_train_x, no_clusters, algorithm="Hartigan-Wong")
fviz_cluster(cluster_data, data = data_train_x)

Tutaj, aby wybrać anomalie spośród punktów w grupach należałoby dobrać jakąś miarę niepodobieństwa. Określmy przykładową miarę niepodobieństwa jako odległość danej od środka grupy do której została zakwalifikowana.

cluster_centers = cluster_data$centers
data = data_train
data %<>% mutate(group = cluster_data$cluster) %>% mutate(dist = sqrt((V1-cluster_centers[group, 1])^2 + (V2-cluster_centers[group, 2])^2 + (V3-cluster_centers[group, 3])^2 + (V4-cluster_centers[group, 4])^2)) %>% arrange(dist)
data$dist = data$dist / max(data$dist)
plot(data$dist)

Miara niepodobieństwa mapuje nam dane wielowymiarowe na jeden wymiar.

plot_ly(data,
        x = ~V1,
        y = ~V2,
        z = ~V3,
        color=~dist,
        colors=c('#ffff00', '#000000'),
        opacity=0.7,
        size=1
      )

## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter3d

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

Jednowymiarowe dane jesteśmy w stanie porównać, dzięki czemu możemy określić jakiś próg, powyżej którego będziemy dane traktowali jako anomalie. Po zakwalifikowaniu anomalii, możemy porównać wyniki z faktycznie wykrytymi anomaliami podanymi w danych wejściowych. W ten sposób dowiemy się, jaka jest jakość modelu a-priori. Na tym etapie wybieramy jakimi miarami jakości chcielibyśmy się posłużyć. W naszym przypadku wybieramy analizę ROC oraz miarę F1.

outlier_threshold = 0.35
data = data %>% mutate(pred_int = as.integer(dist > outlier_threshold))

pred = prediction(predictions=data$dist, labels=data$class)
perf = performance(pred, "tpr", "fpr")
plot(perf)
abline(0,1,col="#0000ff")

Sprawdźmy jeszcze miarę F1. Zwraca wysokie wyniki tylko, gdy jednoczenie oba współczynniki, precyzji i odzysku, dają wysokie wyniki.

data = data %>% mutate(tp = class & pred_int)
data = data %>% mutate(fp = !class & pred_int)
data = data %>% mutate(tn = !class & !pred_int)
data = data %>% mutate(fn = class & !pred_int)
recall = tally(data, tp) / (tally(data, tp) + tally(data, fn))
precision = tally(data, tp) / (tally(data, tp) + tally(data, fp))
f_value = 2 * recall * precision / (recall + precision)
f_value

Na tym etapie mamy wyniki jakie metody grupowania mogłyby dać dla tego zbioru danych. Są to niestety miary które niewiele mówią o jakości predykcji dla nowych danych. Są one w zasadzie tylko poglądowe, aby upewnić się że można wytrenować model w sposób nienadzorowany metodami grupowania, tak aby wyniki detekcji anomalii na podstawie tego modelu były nielosowe i relatywnie poprawne.

W tym miejscu warto zaznaczyć, że wynik NaN oznacza, że i precyzja i odzysk osiągnęły wartość 0, czyli nie było poprawnie zakwalifikowanej nawet jednej anomalii. W pozostałych przypadkach, określona wartość będzie plasowała się w przedziale [0,1], gdzie 0 będzie oznaczało jakość niską, a 1 wysoką.

Natomiast, aby określić jakość modelu, musielibyśmy zastosować go dla nowych danych i sprawdzić czy predykcje okazałyby się identyczne z faktycznym wynikiem kwalifikacji anomalii.

Aby określić jak model radzi sobie z nowymi danymi, będziemy musieli go zastosować dla zbioru danych który nie został wzięty pod uwagę przy tworzeniu modelu.

outlier_threshold = 0.4
predict.kmeans <- function(object, newdata){
    centers <- object$centers
    n_centers <- nrow(centers)
    dist_mat <- as.matrix(dist(rbind(centers, newdata)))
    dist_mat <- dist_mat[-seq(n_centers), seq(n_centers)]
    max.col(-dist_mat)
}

groups = predict(cluster_data, data_test_x)
data = data_test %>% mutate(group = groups)
data %<>% mutate(dist = sqrt((V1-cluster_centers[group, 1])^2 + (V2-cluster_centers[group, 2])^2 + (V3-cluster_centers[group, 3])^2 + (V4-cluster_centers[group, 4])^2)) %>% arrange(dist)
data$dist = data$dist / max(data$dist)
data = data %>% mutate(pred_int = as.integer(dist > outlier_threshold))
plot(data$dist)

Zrobiliśmy predykcje grup na podstawie modelu, zaaplikowaliśmy miarę niepodobieństwa na wyniki danych testowych. Rozkład posortowanych miar niepodobieństwa dla danych testowych prezentuje się jak wyżej.

Aby określić jakość predykcji, zastosujemy analizę roc oraz wyliczymy miarę F1.

pred = prediction(predictions=data$dist, labels=data$class)
perf = performance(pred, "tpr", "fpr")
plot(perf)
abline(0,1,col="#0000ff")

data = data %>% mutate(tp = class & pred_int)
data = data %>% mutate(fp = !class & pred_int)
data = data %>% mutate(tn = !class & !pred_int)
data = data %>% mutate(fn = class & !pred_int)
recall = tally(data, tp) / (tally(data, tp) + tally(data, fn))
precision = tally(data, tp) / (tally(data, tp) + tally(data, fp))
f_value = 2 * recall * precision / (recall + precision)
f_value

Miary jakości powyżej oznaczają, że detekcja anomalii z użyciem takich parametrów, miary niepodobieństwa w postaci odległości od najbliższego klastra i z zastosowaniem grupowania centroidami nie przynosi znaczącego uzysku ponad losowe wybieranie anomalii spośród danych. Może być to spowodowane, że algorytm centroidów jest słabo dostosowany do charakterystyki danych tego zbioru, jak również słabym dopasowaniem miary niepodobieństwa.

Możliwe, że w przypadku tego zbioru danych zastosowanie innych algorytmów grupowania oraz innych miar niepodobieństwa do grup okazałoby się jakościowo lepszym rozwiązaniem.

Klasyfikacja

Do określania miary F1 użyjemy definicji funkcji:

f_value <- function (cm) {
  fp = cm[1,2]
  fn = cm[2,1]
  tp = cm[2,2]
  recall = tp / (tp+fn)
  precision = tp / (tp+fp)
  f = (2 * recall * precision) / (recall + precision)
  return(f)
}

Algorytmy

Support Vector Machines (SVM)

svm.model = svm(as.factor(class) ~ ., data = data_train, type = 'C-classification', kernel = 'polynomial', probability = TRUE)
svm.pred = predict(svm.model, newdata = data_test_x)
svm.cm = table(data_test_y$class, svm.pred)
f_value(svm.cm)

## [1] 0.6086957

Naive-Bayes

nvb.model = naiveBayes(as.factor(class) ~ ., data = data_train)
nvb.pred = predict(nvb.model, newdata = data_test_x)
nvb.cm = table(data_test_y$class, nvb.pred)
f_value(nvb.cm)

## [1] 0.3700787

C-Tree

ct.model = ctree(as.factor(class) ~ ., data = data_train)
ct.pred = predict(ct.model, newdata = data_test_x)
ct.cm = table(data_test_y$class, ct.pred)
f_value(ct.cm)

## [1] 0.5511811

Jakość modeli

require(ROCR)
svm.prob = predict(svm.model, type="prob", newdata= data_test_x, probability = TRUE)
svm.rocr = prediction(attr(svm.prob, "probabilities")[,2], data_test_y$class)
svm.perf = performance(svm.rocr, "tpr","fpr")
nvb.prob = predict(nvb.model, newdata=data_test_x, type="raw")
nvb.rocr = prediction(nvb.prob[,2], data_test_y$class)
nvb.perf = performance(nvb.rocr, "tpr","fpr")
ct.prob = predict(ct.model, newdata=data_test_x, type="prob")
ct.rocr = prediction(ct.prob[,2], data_test_y$class)
ct.perf = performance(ct.rocr, "tpr","fpr")
plot(svm.perf, col=4, main="ROC curves of different machine learning classifier")
legend(0.6, 0.6, c('svm', 'ctree', 'Naive Bayes'), 4:6)
plot(ct.perf, col=5, add=TRUE)
plot(nvb.perf, col=6, add=TRUE)